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In material sciences, spectroscopic approaches combining ab initio calculations with experiments 
are commonly used to accurately analyze the experimental spectral data. Most state-of-the-art 
first-principle calculations are usually performed assuming an equilibrium static lattice. Yet, nuclear 
motion affects spectra even when reduced to the zero-point motion at 0 K. We propose a framework 
based on Density-Functional Theory that includes quantum thermal fluctuations in theoretical X- 
ray Absorption Near-Edge Structure (XANES) and solid-state Nuclear Magnetic Resonance (NMR) 
spectroscopies and allows to well describe temperature effects observed experimentally. Within the 
Born-Oppenheimer and quasi-harmonic approximations, we incorporate the nuclear motion by gen¬ 
erating several non-equilibrium configurations from the dynamical matrix. The averaged calculated 
XANES and NMR spectral data have been compared to experiments in MgO, proof-of-principle 
compound. The good agreement obtained between experiments and calculations validates the de¬ 
veloped approach, which suggests that calculating the XANES spectra at finite temperature by 
averaging individual non-equilibrium configurations is a suitable approximation. This study high¬ 
lights the relevance of phonon renormalization and the relative contributions of thermal expansion 
and nuclear dynamics on NMR and XANES spectra on a wide range of temperatures. 


I. INTRODUCTION 

X-ray Absorption Near-Edge Structure (XANES) ^ and 
solid-state Nuclear Magnetic Resonance (NMR)^ spec¬ 
troscopies are powerful probes of the electronic and lo¬ 
cal structure of inorganic materials. The combination 
of these two techniques provides a deep understanding 
of the electronic and structural properties of materials. 
For instance, a recent study coupling X-ray absorption 
spectroscopy and NMR successfully resolved the local 
structure of A1 sites in zeolites.^ The noticeable improve¬ 
ments in methodology and instrumentation resulted in 
a spectacular enhancement of the quality and resolution 
of spectra for both techniques.Nonetheless, the huge 
amount of information contained in experimental spec¬ 
tra makes their accurate assignment difficult. To address 
this difficulty theoretical tools are used beforehand or in 
conjunction with experimental data.®“^° 

Most first-principle calculations in the solid-state con¬ 
sider the nuclei fixed at their equilibrium positions, as ob¬ 
tained by X-ray or neutron diffraction. However, atoms 
are subjected to quantum thermal fluctuations, which 
reduce to the zero-point motion at the absolute zero 
(0 K). Temperature-dependent experiments exhibit sig¬ 
nificant variations of spectroscopic properties, such as 
chemical shifts,relaxation times,and Electric Field 


Gradient values^^^^^ (EFG) measured by NMR or the 
pre-edge structures observed in XANES spectra of oxide 
materials.In the case of NMR, it has been observed 
that the chemical shifts vary of several ppm when tem¬ 
perature increases over thousands of degrees.xhe 
EFG tensor is also very sensitive to temperature and 
nuclear quadrupole resonance experiments have shown 
a variation of the resonances frequencies o up to 
—0.2 kHz.K“^.^^d5 jji case of XANES spectroscopy, 
it has been recently shown that the intensity and position 
of the pre-edge peak highly depend on temperature. 

It has been demonstrated that the pre-edge structure is 
due to a violation of symmetry induced by the quantum 
thermal fluctuations. For instance, at the A1 A'-edge 
of various oxides, the nuclear motion is responsible for 
the appearance of a pre-edge peak corresponding to the 
Is —>■ 3s forbidden electronic transitions. 

Multiple attempts to theoretically reproduce the lat¬ 
tice dynamical effects in solid-state spectroscopies have 
emerged in the literature. Pioneering approaches based 
on averaging the chemical shielding tensors over differ¬ 
ent orientations of mobile species were proposed to in¬ 
clude nuclear motion in NMR calculations.In par¬ 
allel, it was shown that small displacements of either 
the absorbing atom^°’^^’^® or the Is initial wave func¬ 
tion in the crude Born-Oppenheimer approximation^^ 
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could reveal forbidden transitions in K-edge XANES 
spectra. Although quite promising, these methods do 
not account for the collective lattice dynamics. A sub¬ 
stantial theoretical work proved that vibrations could 
be represented as the convolution of the x-ray absorp¬ 
tion cross section, calculated for the equilibrium con¬ 
figuration, with the phonon spectral function.This 
theory has been applied to reproduce room-temperature 
experimental XANES spectra in Ref. 20. In the case 
of non-quantum nuclear motion, the method of choice 
would be Molecular Dynamics (MD) at finite temper¬ 
ature, either classical or ab initio. In NMR, ab initio 
MD was used to study the dependence on temperature 
of the chemical shift^^’^®“^^ and quadrupolar relaxation 
rates^^’^^ but mostly in organic compounds. MD cal¬ 
culations were also used in XANES to calculate Li and 
S AT-edges XANES spectra of Li-ion batteries at room 
temperature^^’^^ or the A1 and Fe AT-edges in warm dense 
plasmas.However, MD consider the vibrations as a 
classical phenomenon and is therefore appropriate only 
if ksT > hzuyib, with zuvib the vibration frequency.^® To 
account for the quantum behavior of vibrations. Path 
Integral Molecular Dynamics (PIMD) was used to simu¬ 
late NMR^® and C AT-edge NEXAFS^*^ spectra in organic 
compounds. In both spectroscopies, the results were im¬ 
proved using PIMD, but at the cost of a larger statis¬ 
tical error. Seeking for computationally less expensive 
methods, Monte Carlo sampling has been used to ac¬ 
count for the vibrational effects on the chemical shifts 
in MgO^^ and recently a more computationally efficient 
method arose by assuming a quadratic coupling between 
vibrations and the shielding tensor.The Monte Carlo 
sampling method was also used to simulate XANES spec¬ 
tra of molecules in solutions,but, to our knowledge, it 
has not been applied to inorganic solids yet. 

The purpose of this work is to describe quantum 
thermal fluctuations, using methods based on Density- 
Functional Theory (DFT).^^ In the Born-Oppenheimer 
(BO)^® and Quasi-Harmonic Approximation (QHA),^®d7 
the thermal effects are modeled by generating atomic 
position configurations obeying quantum statistics at fi¬ 
nite temperatures. The theoretical work is confronted 
to NMR and recently acquired XANES temperature- 
dependent experiments for our proof-of-principle com¬ 
pound MgO. This ionic oxide has been chosen for three 
main reasons: (i) its rock salt structure is ideal to ob¬ 
serve thermal effects as its lattice constant is the only 
free parameter, (ii) it shows no phase transition up to its 
melting point around 3250 K,^® allowing a wide tempera¬ 
ture range for experimental measurements, (iii) multiple 
theoretical studies demonstrated that the harmonic be¬ 
havior of MgO remains at temperature as high as 1500 K 
and its Debye temperature has been found, theoretically, 
to be about 941 K.^s-si 

The paper is structured as follows. In Sec. H, the for¬ 
malism upon which our theoretical model is based is de¬ 
tailed. In Sec. HI, experimental and computational de¬ 
tails are given. In Sec. IV, experimental and theoretical 


results obtained on MgO are presented and discussed. 
Finally, the conclusions of this work are drawn in Sec. V. 

II. FORMALISM 

A. Quasi-harmonic vibrations 

A crystal can be seen as a system of N nuclei and 
electrons with respective position vectors (Ri,..., Rat) 
and (ri,..., ttv,,). The collective coordinates R = 
(Ri,..., Rat) and r = (ri,..., tat,,) are used thereafter. 
The stationary states of the system are described by the 
wave function 4'(r, R) solution of the following general 
Schrddinger equation, 

(Tat + Te + Ve + Ve-N + Vat) ^(r,R) = A^(r,R), ( 1 ) 

with Tat the kinetic nuclear operator, Tg the kinetic elec¬ 
tronic operator, 14 the Coulomb potential between elec¬ 
trons, Ve-N the Coulomb potential between nuclei and 
electrons, and I4v the Coulomb potential between nuclei. 
The total energy of the crystal is denoted by E. 

In the BO approximation, which assumes that the elec¬ 
tronic cloud reacts instantaneously to the nuclear motion, 
the wave function solution of Eq. (1) can be approxi¬ 
mated as the product 

4/^„(r,R)=x^„(R)V^„(r;R), (2) 

of an electronic part R), in which R is a parameter, 
and a nuclear part (R), xhe electron orbital index is 
n, and j indexes vibrational states. 

The Hamiltonian acting onto the electronic variables 
is labelled Hbo, 

HbO = Te + Ve + Ve-N + Vn ( 3 ) 

where Vn is then a constant energy term determined for a 
given nuclei configuration R. Thus Hbo is parametrized 
by R. The electronic wave function '0„(r;R) is solution 
of the Schrddinger equation, 

HBOipniv R) = e„(R)4>„(r; R), (4) 

which introduces the energy surface £:„(R). In the BO 
approximation the lattice dynamics is described by 

[Tn + £„(R)] xi(R) = EixiW (5) 

where the vibrational wave functions (R) are the or¬ 
thonormal solutions and is the total energy of the 
whole crystal. Equations (4) and (5) imply that {ipn} 
and {xn} are complete basis sets of eigenvectors of Hbo 
and of the nuclear Hamiltonian [Tat-|- en(R)], respec¬ 
tively, leading to the completeness relations 

=S(r-r');'iR, (6) 

n 

X^„(R)x^n(R') = S(R - R'); Vn. (7) 
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The phonon-induced displacement of nucleus / in the 
Cartesian direction a 

u^=RJ- (8) 

is small compared to the atomic bond length. Thus a 
Taylor expansion of the BO energy surface as a function 
of the nuclear displacements can be carried out and is 
truncated at the second order in the harmonic approxi¬ 
mation. The energy scale is shifted such as the zero-order 
term is null and the first-order term vanishes because 
forces acting on individual nuclei are zero at equilibrium. 
In this approximation, [Tjv-I-en(R)] becomes the nu¬ 
clear harmonic Hamiltonian as 
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as a product of Gaussian functions following a normal 
distribution whose widths depend on temperature and 
phonon frequency.^^ The Gaussian functions are centered 
on = 0, i.e., at the equilibrium position when u/ = 0 
(Eq. 12). The 7^(R) probability distribution is written 
as 


P(R) = ylexp 
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with A a normalization constant. The normal length 
of the ^ vibration mode is the standard deviation of col¬ 
lective normal coordinates q and depends explicitly on 
temperature T and phonon frequencies 



(16) 


where Tv is explictly written in terms of P/, the mo¬ 
mentum operator of the /th nucleus. In Eq. (9) we have 
introduced the interatomic force constant matrix C whose 
elements are 


qIJ^ _ 


9^en(R) 


dufduj 


( 10 ) 


(eq) 


Rescaling C by the nuclei masses define the dynamical 
matrix, whose diagonalization as 


N 3 


EE 

J—l — l 




( 11 ) 


provides phonon polarization vectors e“^ and phonon fre¬ 
quencies where indexes phonon modes. Using in 
Eq. (9) 
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which introduce the normal coordinates and = q^.-, 
% can be written as a sum of 3iV independent Hamilto¬ 
nian operators of harmonic oscillators 


with P = l/ksT. Equation (16) describes how the 
phonon normal modes are thermally populated. The sta¬ 
tistical average of any observable O is then 

{O) = J dR 0(R) T(R)- (17) 


For a given crystal, after calculating and diagonaliz¬ 
ing the dynamical matrix (Eq. 11) to obtain the phonon 
frequencies and polarization vectors e^, a set of Nc 
nuclear configurations obeying the T(R) quantum sta¬ 
tistical distribution is generated. For each p mode, a 
set of random Gaussian numbers is cre¬ 

ated. Then, each is multiplied by the corresponding 
normal length a^, as defined in Eq. (16). The set of so- 
generated normal coordinates {g), = obeys 

the probability distribution (Eq. 15). The nuclear po¬ 
sition collective vectors {r2-i.. are obtained using 
Eqs. (8) and (12) such as 




= R 
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According to the importance sampling technique, 
Eq. (17) is equivalent to 


(O) 


1 






(19) 
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77 = ^ 2 ^ 17 ) ■ ( 14 ) 
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Nevertheless, a harmonic model considers the phonon 
normal modes as independent quasiparticles and does not 
account for any anharmonic phenomenon, such as ther¬ 
mal expansion. In this work QHA is used to include 
thermal expansion effects. The model is no longer purely 
harmonic but does not describe phonon-phonon interac¬ 
tion, as the phonon normal modes are still independent. 
Within QHA, the probability T(g^) of finding the sys¬ 
tem in any set of normal coordinates g^, is expressed 


In this work, 0(R*) is either the nuclear magnetic shield¬ 
ing tensor, the EFG tensor or the XANES cross section 
for the Ah configuration. 

B. Nuclear Magnetic Resonance spectroscopy 

To compare with NMR experiments, for each nucleus, 
the isotropic value of the magnetic shielding tensor cr is 
evaluated as 

1 / ^ 

CTiso = 3 l''('^) • 


( 20 ) 
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More precisely, the isotropic chemical shift Siso — 
— (ciso ~ CTref) is considered, where a^ef is the isotropic 
shielding of the standard reference. Unlike the Larmor 
frequency, (5iso contains the response of the electronic sys¬ 
tem to the magnetic perturbation and is therefore an ex¬ 
plicit function of the electronic environment around the 
probed nucleus. Thus, 5iso = <5iso(R-) and ((5iso(R.)) is 
obtained as in Eq. (19). 

C. XANES spectroscopy 

In a single-electron approach, at the K-edge, the 
XANES cross-section, in the electric dipole approxima¬ 
tion is usually given by^^ 

cr(ftw) = ® - eis - fiw). 

n 

( 21 ) 

where fej and e are the energy and polarization vector 
of the incident X-ray photon, respectively, and oq the 
fine structure constant. The final \ipn) and initial \ijjis) 
electronic states of energy £„ and eis, respectively, are 
solution of Eq. (4). 

Equation (21) rests upon the assumption that the 
nuclei are fixed in a generic configuration, hence a is 
parametrized by R, i.e., aihjj) = (7(/ia;;R). A more rig¬ 
orous expression of the cross section accounting for the 
total system of nuclei and electrons is 

crtot(Sw) = dTT^ao^ia;^ - f 

n.j 

( 22 ) 

where the initial state and final states are defined 
as in Eq. (2). In Appendix A, it is demonstrated that atot 
can be expressed, using the first order of the expansion 
given in Eq. (A3), by 

crtot{f^) = y dR p(R)cr(fiw;R), (23) 

with p(R) the quasi-harmonic weighting displacement 
distribution. The next orders of the expansion in 
Eq. (A3) are not considered in the present work. In 
Eq. (23), atot is the average of individual configuration 
cross sections a{fiuj; R) using a probability distribution 
that takes into account the temperature and vibrations 
frequencies in a form consistent with Eq. (19). 

III. EXPERIMENTAL AND CALCULATION 
DETAILS 

A. Experimental setup 

Mg AT-edge X-ray absorption experiments were per¬ 
formed at the LUCIA beam line (SOLEIL Saint-Aubin, 
France).®^ The incident energy range (1280-1400 eV) was 
selected to include the Mg A-edge using a double Beryl 


monochromator. The pressure in the experimental cham¬ 
ber was 10“^ mbar. The incident X-ray beam was set to 
a spot size of I x 2 mm^. Temperature-dependent mea¬ 
surement were conducted at 300 K, 573 K, 773 K, 873 K 
and 1273 K using a boron nitride furnace. Only the spec¬ 
tra recorded at 300 K, 573 K and 873 K are shown. The 
4 cm^ MgO single crystal was held using a perforated 
lamella of molybdenum. The temperature of the sam¬ 
ple was measured using a Chromel-Alumel thermocouple. 
The spectra were recorded in fluorescence mode with a 
four element Silicon Drift Diode detector, protected from 
infrared and visible radiations by a thin beryllium win¬ 
dow. To maximize the signal/noise ratio, each point was 
obtained after a 5 second acquisition time and 5 spectra 
were measured for each temperature. The self-absorption 
correction and spectra normalization were applied as in 
Ref. 20. 

The temperature-dependent measurements of ^®Mg 
and static isotropic chemical shifts in MgO are taken 
from Ref. 11. The NMR active isotopes ^®Mg and 
have weak natural abundances and gyromagnetic ratios, 
thus experiments required isotopically enriched samples. 
Both nuclei are quadrupolar (/spin = 5/2) but, as any 
atomic site in MgO presents Oh symmetry, the exper¬ 
imental EFG vanishes. Therefore, NMR peaks do not 
suffer from quadrupolar broadening and shifting. 


B. Calculation details 


TABLE I. Parameters of generation for Troullier-Martin ul- 
trasoft (US) and norm-conserving (NC) pseudopotentials. 
Bessel functions are used to pseudize the augmentation 
charges. The radii are in Bohr units. 


Atom 

Valence states (Radius) 

Local part 

Mg (NC) 

3s^(2.00) 3p°(2.00) 3d°(2.00) 

d 

O (NC) 

2s^(1.45) 2p®(1.45) 

p 

Mg" (US) 38^(2.50) 3p°(2.60) 3d®(2.30) 

d 

0 (US) 

2sU1.35) 2p'^(1.35) 

p 


The pseudopotential of the absorbing Mg atom was generated 
using the same parameters but with only one Is core electron 
and used in the XANES calculation. 


All the calculations were performed using the pseu¬ 
dopotential, plane wave Quantum ESPRESSO suite of 
codes,®® within the DFT-PBE generalized gradient ap¬ 
proximation (GGA).®® The details of the pseudopoten¬ 
tial used herein are given in Table I. Most of the calcu¬ 
lations were done using ultrasoft®^ GIPAW®® pseudopo¬ 
tentials except for the NMR calculations for which norm- 
conserving Troullier-Martin®® GIPAW pseudopotentials 
were preferred. 

The method detailed in Sec. II has been carried out 
creating the configurations at the Quasi-Harmonic level 
with the Stochastic Self-Gonsistent Harmonic Approxi- 
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mation (SSCHA) code®°’®^ starting from the QHA dy¬ 
namical matrices of MgO, which crystallizes in space 
group Fm3m with a room temperature lattice parame¬ 
ter a = 4.21 A. The NMR and XANES calculations were 
performed for temperatures ranging from 0 K to 1273 K 
with the experimental lattice parameters,®^ so that the 
configurations and spectroscopic results account for the 
thermal-expansion anharmonic effect. The 0 K calcula¬ 
tions were performed using the 12 K lattice parameter.®® 
Self-consistent electronic densities were calculated in 
the unit-cell at the volume of each corresponding tem¬ 
perature. A 4 X 4 X 4 k-point grid sampled the Bril- 
louin zone,®^ and a plane-wave energy cutoff (density) 
for the wave functions of 45 (540) Ry was chosen. Then, 
the dynamical matrices were calculated using Density- 
Functional Perturbation Theory®®’®® on a q-point grid 
commensurable with the supercell size needed for NMR 
and XANES calculations thereafter. Long-range electro¬ 
static interactions were taken into account by calculating 
Born effective charges and electronic dielectric tensor.®^ 
NMR calculations were performed using the GIPAW 
package.®®’®® Eor each temperature, the calculations were 
conducted for 2x2x2 supercell configurations containing 
64 atoms. Each calculation was done on a 2 x 2 x 2 k- 
point grid, with a plane-wave energy cutoff for the wave 
functions (resp. density) of 90 (resp. 360) Ry. The con¬ 
vergence was reached for ten configurations at each tem¬ 
perature. Eor a given temperature, 640 isotropic shield¬ 
ing tensors were calculated and averaged for each nu¬ 
cleus. In this work, Uref is chosen so that experimental 
and calculated values match at room temperature. The 
principal components Vxx, Vyy and Vzz of the EEG ten¬ 
sor defined with IVzzj > |I4x| > |Vyy| are obtained by 
diagonalization of the tensor. The quadrupolar coupling 
constants Cq defined as Cq = eQVzz/h were estimated 
with Q values from Ref. 69 and were found negligible at 
each temperature. This behavior is consistent with the 
crystal symmetry and experimental results. Indeed, in a 
perfect cubic environment we expect the first-order EEG 
tensor quantities to vanish. 

Eor each temperature, the XANES spectra were calcu¬ 
lated using the XSpectrapackagein 3x3x3 supercell 
configurations containing 216 atoms, each including one 
Is full core-hole in a random Mg site. The self-consistent 
electronic density was obtained at the T point of the Bril- 
louin zone with a plane-wave energy (resp. density) cut¬ 
off of 45 (resp. 540) Ry. The XANES theoretical spectra 
were performed on a 4 x 4 x 4 k-point grid with a con¬ 
stant broadening parameter of 0.5 eV. Eor a given con¬ 
figuration, three polarized XANES spectra were calcu¬ 
lated for the X-ray polarization vector e parallel to each 
of the Gartesian axes and the average gave an isotropic 
XANES spectrum. As a convergence criterium, it has 
been verified that the averaged polarized spectra along 
each Cartesian direction matched, as expected in a cubic 
crystal. The convergence was reached for a number of 
30 configurations at each temperature. In the pseudopo¬ 
tential approximation, only valence-electrons are consid¬ 


ered, hence the energy scale has no physical meaning and 
a core-level shift has to be applied before comparing and 
averaging theoretical spectra (Eq. 19). Similarly to pre¬ 
vious works,the core-level shift in the ith configu¬ 
ration is taken into account as follows: 

. (24) 

In this rescaling, the energy of the lowest unoccupied 
electronic band was subtracted and the energy dif¬ 
ference between the system with one Is core-hole and one 
electron in the first available electronic state and 

that of the ground state (E^g) was added. Einally, all 
the spectra were shifted by 1303.8 eV to match the ex¬ 
perimental main edge peak energy position at room tem¬ 
perature. To interpret XANES spectra, local and partial 
density of states (DOS) were calculated in the 3x3x3 
supercell configurations using Lowdin projections on a 
4x4x4 k-point grid with a Gaussian broadening pa¬ 
rameter of 0.3 eV. 


IV. RESULTS AND DISCUSSION 
A. Solid-state NMR spectroscopy 

Eigure 1 compares the temperature-dependence of the 
experimental NMR measurements from Ref. 11 with the 
calculated isotropic chemical shifts i5iso of ^®Mg and 
in MgO up to 1273 K. Eor both nuclei, the calculated val¬ 
ues closely reproduce the experimental trend: the chem¬ 
ical shift increases with temperature. In addition, the 
contribution of the thermal expansion is displayed. The 
incorporation of the quantum motion of nuclei is manda¬ 
tory, since it improves the agreement with experiment, 
especially for ^®Mg. Without the quantum motion, con¬ 
sidering only the thermal expansion leads to a downward 
trend in ^®Mg and an upward one in These opposite 
trends come from the fact that unlike ^®Mg, the chemical 
shift of is not a function of the Mg-0 bond length. 
Indeed, the chemical shift of depends on the inter¬ 
action between the empty d states of ^®Mg and the 2p 
states of which results in a deshielding of ®^0.^® 

TABLE II. Slopes of the temperature-dependent chemical 
shifts (5iso for temperatures ranging from 300 K to 1300 K, 
in ppm.K~^. The given values are the best fit lines in the 
room-to-high temperature region. The experimental results 
are reported and compared to calculated values from this work 
and previous studies. 


Nucleus 

Expt.^ 

This work 

Calc.® 

Calc." 

2®Mg 

0.002 

0.003 

0.004 

0.005 


0.008 

0.010 

0.011 

0.005 


Ref. 11 
® Ref. 41 
Ref. 42 
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FIG. 1. (Color online) Temperature-dependent isotropic chemical shift of ^®Mg and in MgO obtained experimentally^^ 
(diamonds), calculated in the QHA (circles), and considering only the thermal expansion (triangles). The linear fits of each 
dataset are displayed. For the QHA calculation, two independent linear contributions are observed: from 0 K to 100 K and 
from 300 K to 1273 K. The experimental errors bars are set to ±1 ppm as in the original paper, while the calculated error bars 
are the statistical uncertainties. 


The linear trend behavior of the chemical shift ob¬ 
served from 300 K to 1273 K is consistent with exper¬ 
iments and previous theoretical studies.Our results 
show a better agreement with experiments than other 
theoretical studies, the corresponding slopes obtained for 
each nucleus are summarized in Table II. The method 
used in Ref. 41 is similar to ours, but the remaining dis¬ 
crepancies may come from a less accurate description of 
the phonons dispersion and quasi-harmonic vibrational 
wave functions. In Ref. 42, the coupling between the 
phonons and the chemical shift tensor was expanded in 
terms of vibrational-mode amplitudes and was assumed 
to be quadratic. However, in this study, the thermal ex¬ 
pansion was neglected, thus, the results cannot be com¬ 
pared to experiment and the same slope is reported for 
the temperature-dependence of both nuclei unlike exper¬ 
iments. 

The linear behavior in the room-to-high temperature 
region arises from a combination of both the thermal 
expansion and dynamics of nuclei. These two effects 
are constant in the low-temperature regime, from 0 K 
to 100 K, where the vanishing thermal expansion gives 
rise to a flat temperature dependence and the dynamics 
of the nuclei reduces to a Ai5iso constant term (Fig. 1). 
The zero-point renormalization of the chemical shift is 
evaluated to A(5iso(^®Mg) = 0.69 ppm and A(5iso(^^0) = 
0.60 ppm, respectively. Above 300 K, the amplitude of 
the dynamical part outweighs the thermal expansion. Fi¬ 
nally, in the intermediate temperature region (typically 
between 100 K and 300 K) there is a competition between 
the two components of the temperature-dependence that 
leads to the observed curvature. 

The origin of the remaining discrepancies between cal¬ 
culated and experimental slopes has to be investigated. 


The overestimation of the temperature-dependence of 
both nuclei chemical shifts could be related to a defi¬ 
ciency of GGA. Indeed, it is well known that GGA over¬ 
estimates slightly the interatomic distances of solids and 
underrates phonon frequencies.^® Although experimental 
lattice parameters are chosen to balance this effect, the 
use of more accurate GGA functionals^^’^® might enhance 
the agreement. Furthermore, it has been noticed that, in 
ionic compounds containing alkaline-earth cations, the 
calculation-experiment correlation lines of NMR param¬ 
eters, as calculated within GGA, deviate from the 1 ideal 
value.In the case of CaO, it has been observed that 
the shortening of the bond length induces an artificial 
hybridization between the 3d states of Ca and the 2p of 
O.^® This problem was addressed by the use of a cor¬ 
rected Ca pseudopotential. To a lesser extent this prob¬ 
lem can occur in MgO. Laskowski et al. proposed to 
use GGA+U instead of standard GGA to better describe 
the 3d states of the cation. Finally, one may argue of 
a failure of QHA, but the theory has revealed to accu¬ 
rately provide the lattice dynamical properties in MgO 
up to 1100 K under ambient pressure conditions.®® To 
go beyond QHA, anharmonic effects in periodic solids 
can be investigated in a numerically feasible manner us¬ 
ing the vibrational self-consistent field method®^ or the 
SSCHA.®! 


B. XANES spectroscopy 

In Fig. 2(a) the experimental Mg K-edge XANES 
spectra are reported. While temperature continuously 
smoothes the XANES features, the P pre-edge peak in¬ 
creases and shifts towards lower energy. To a lesser 
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FIG. 2. (Color online) Mg K-edge XANES experimental (a) 
and QHA-calculated (b) spectra of MgO, for temperature up 
to 1273 K, along with difference of each spectrum with respect 
to the room-temperature one. 


extent, this energy shift is also visible for the A main 
edge peak and the following features. These tempera¬ 
ture effects have been also observed at the A1 itT-edge in 
corundum (a-Al 203 ) and beryl (Be3Al2Si60i8).^° The 
corresponding calculated XANES spectra are plotted in 
Fig. 2(b) along with the 0 K and 1273 K spectra. Cal¬ 
culations reproduce closely experiments over all the ex¬ 
plored incident X-ray photon energy range, as highlighted 
by the similarity of the difference of each spectrum with 
the room-temperature one. However, the intensity of the 
calculated pre-edge is underestimated: the experimental 
and calculated P/A intensity ratios are about 1/5 and 
1/7, respectively. This mismatch may come from the 
first-order expansion of the X-ray absorption cross sec¬ 
tion (Eq. 23). Going further than the first-order is a 
possible improvement of the method. In addition, as for 
NMR, conhgurations obtained beyond QHA could lead 



FIG. 3. (Color online) Experimental and calculated energy 
positions of P (upper panel) and A (lower panel) peaks as a 
function of temperature. The energy positions are given with 
respect to the A peak position at 300 K. 


to a better agreement with experiment. 

The temperature-dependence of the P and A peaks en¬ 
ergy positions is showed in Fig. 3. For temperature rang¬ 
ing from 300 K to 1273 K, the theory-experiment agree¬ 
ment is satisfactory: the experimental and calculated P 
peak variation is 0.55 eV and 0.41 eV, respectively. For 
peak A, the experimental (resp. calculated) variation is 
0.17 eV (resp. 0.12 eV). The results are comforted by the 
experimental observations at the A1 AT-edge in corundum 
where the pre-edge position decreased of about 0.4 eV 
from 300 K to 930 The origin of these shifts may 
be related to the band gap evolution in temperature. In 
ionic compounds, the lattice expansion affects electronic 
bands. In MgO, the band gap was shown to decrease by 
0.91 eV, from 300 K to 1273 K, using optical reflectivity 
measurements.®^ Moreover, the electron-phonon interac¬ 
tion contributes more than the thermal expansion to the 
band gap narrowing.®®’®^ Over the same range of tem¬ 
perature, our calculations achieved a similar band gap 
decrease (0.78 eV), while the decrease only due to ther¬ 
mal expansion is 0.3 eV. 

The energy shift to lower energy does not exclusively 
concern the P and A features since it is visible over all 
the spectral energy range. For instance, peak D moves of 
about 0.1 eV to lower energy as temperature increases. 
The decreasing shift agrees with the empirical Natoli’s 
rule {E(P = cste),®® which states that the energy position 
decreases with increasing interatomic distance. There¬ 
fore, this signature can be related to the thermal expan¬ 
sion. Figure 4 displays the Mg AT-edge XANES spectra 


























FIG. 4. (Color online) Calculated Mg IF-edge XANES spectra 
in MgO considering only the thermal expansion. 


calculated in the equilibrium configuration at the vol¬ 
umes corresponding to 300 K and 1273 K. The 1273 K 
spectrum is more contracted than the 300 K spectrum. 
The C and D features move down toward lower energy 
with increasing interatomic distance as observed exper¬ 
imentally. On the contrary, in opposite trend with ex¬ 
periment, peak A shifts toward higher energy. There¬ 
fore, thermal expansion does not fully explain the shift¬ 
ing trends, especially at lower energy, where vibrations 
are mandatory to reproduce the correct spectral-feature 
positions. 

Former studies already proposed that the P peak orig¬ 
inates from the vibration-induced violation of the dipole 
forbidden Is —>■ 3s transitions.To further ana¬ 
lyze the P peak origin, Fig. 5 compares the theoretical 
XANES spectra obtained at 0 K, by including or not the 
zero-point motion. In addition. Fig. 5 displays all the 
core-level shifted individual configuration spectra in the 
background. The phonon influence on the XANES spec¬ 
tra is characterized by two main features. First, similarly 
to a convolution and in agreement with the theoretical 
framework of Fujikawa,^® the inclusion of the 0 K quan¬ 
tum fluctuations globally smoothes the XANES spec¬ 
trum. Second, the zero-point motion induces a pre-edge 
peak that is totally absent in the equilibrium spectrum, 
i.e., when the atoms are fixed at their equilibrium posi¬ 
tions. Hence, Fig. 5 highlights the quantum origin of the 
pre-edge. Moreover, the weak difference in the pre-edge 
intensities between 0 K and 300 K calculations empha¬ 
sizes the prominent role of quantum effects up to room 
temperature [Fig. 2(b)]. A better description of the P 
and A peaks intensities and variations could be achieved 
by the phonon-renormalization of the electronic energies 
directly in the self-consistent calculation, as performed 
in Refs. 86-89. 

The 0 K individual configuration spectra present a 
strong dispersion around the averaged spectrum (Fig. 5). 
This dispersion increases with temperature. Some of the 
individual configuration spectra exhibit a strong P peak 
while some other do not, whatever the temperature. Two 



FIG. 5. (Color online) Impact of the zero-point motion on 
the theoretical Mg A-edge XANES spectrum of MgO. The 
dashed spectrum is obtained with the atoms fixed at their 
equilibrium positions in the 12 K experimental volume.The 
solid spectrum is the average of the 30 configuration spectra 
(in light blue). 


individual configuration XANES spectra at T = 1273 K, 
one with a pre-edge and one without, have been selected 
and analyzed using local and partial DOS [Figs. 6(b,c)]. 
For comparison. Fig. 6(a) displays the case of the stan¬ 
dard calculation with atoms fixed at their equilibrium 
positions at T = 1273 K. The partial and local DOS 
plotted in Fig. 6 are the s and p empty DOS projected 
on (i) the absorbing Mg, (ii) the first O neighbors, (iii) 
the next Mg neighbors, and (iv) the next O neighbors. 
Whatever the configuration, the core-hole strongly mod¬ 
ifies the s and p empty states of Mg, leading to two 
rather localized peaks, which coincide with P and A, re¬ 
spectively. P is only visible if sp hybridization of the 
Mg absorbing states occurs [Figs. 6(c)]. This sp mix¬ 
ing of the Mg absorbing states also induces a stronger 
p DOS of oxygen. The sp hybridization due to the dy¬ 
namical distortion of the MgOe octahedron does not sys¬ 
tematically occur [Fig. 6(b)]. The contribution of the 
lowest unoccupied band to the electronic charge density 
(± I'i/'lub]^) gives a representation of the electronic state 
probed in the P peak and enlightens the effect of the 
MgOe octahedron distortion. In the equilibrium configu¬ 
ration [Fig. 6(a)] the isosurface shows a centrosymmetric 
cubic shape. The distortion of the lattice strongly im¬ 
pacts ± IV’lub] [Fig- 6(b)], however, it is not sufficient to 
create a pre-edge peak. The P peak emerges if the distor¬ 
tion of the MgOe octahedron induces a p-like character 
on the neighboring O atoms, as already observed in the 
DOS [Fig. 6(c)]. To conclude, the breakdown of the sym¬ 
metry is mandatory to induce the sp hybridization and 
the forbidden Is —> 3s transition but is not sufficient. 
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FIG. 6. (Color online) Pre-edge peak analysis. Three configurations at 1273 K are considered: the equilibrium structure (a), 
an individual configuration not leading to pre-edge peak (b) and an individual configuration that leads to pre-edge peak (c). 
Top panel: Mg iF-edge theoretical XANES spectra. Middle panel: local and partial electronic density of states (see text for 
details). Bottom panel: isosurfaces of sign(r/)LUB) I'J/'lubI^, where r/iLue is the electronic wave function of the lowest unoccupied 
band, whose energy is indicated by the vertical dashed line in top and middle panels (negative and positive signs are colored 
in green and blue, respectively). The Mg atoms are displayed in orange and O atoms in red. The pictures are centered on the 
absorbing Mg atom. The isosurface level is set to 6 x 10“"^ with ao the Bohr radius. 


V. CONCLUSION 

A DFT-based approach enabling to successfully in¬ 
troduce quasi-harmonic quantum thermal fluctuations of 
nuclei in NMR and XANES spectroscopies has been pre¬ 
sented. This method, avoiding the explicit calculation 
of the electron-phonon coupling parameters, provides an 


efficient framework to analyze phonon effects occurring 
in both spectroscopies. The calculated spectral data ob¬ 
tained in the MgO proof-of-principle compound are in 
good agreement with experimental datasets, which sup¬ 
ports the reliability of our approach. 

The combination of experiments and first-principle cal¬ 
culations have enabled to investigate the influence of the 
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quantum vibrations in both spectroscopies. A similar be¬ 
havior is revealed: the zero-point phonon renormalization 
of NMR and XANES spectra improves the experiment- 
calculation agreement and therefore could be used on a 
regular basis to reproduce experimental data even at low 
temperatures. 

In NMR, the experiment-calculation agreement is im¬ 
proved with respect to previous theoretical studies. The 
temperature-dependence of the chemical shifts results 
from both contributions of thermal expansion and nu¬ 
clear dynamics, and reduces to a constant renormaliza¬ 
tion term at low temperature. 

In the case of XANES, the temperature-dependence 
of XANES features is reproduced over a large range of 
temperatures. An analytic expression of the phonon- 
dependent X-ray absorption cross section is derived at 
the first order of the electronic Green’s function expan¬ 
sion. It appears that keeping only the first term in the 
expansion of Eq. (A3) is a suitable approximation to 
calculate XANES spectra at finite temperature. It ap¬ 
pears that the first-order calculation of XANES spectra 
at finite temperature is a suitable approximation. Nev¬ 
ertheless, the implementation of the higher-orders cor¬ 
rection terms could improve the pre-edge intensity mod¬ 
eling. The presence of the pre-edge feature is a rele¬ 
vant signature of phonon effects. A thorough study of 
the mechanism from which the pre-edge emerges is con¬ 
ducted. The breakdown of the coordination symmetry is 
mandatory to induce the pre-edge and a p-like character 
arises on the neighboring O atoms if the pre-edge feature 
is discernible. The pre-edge energy variation in tempera¬ 
ture is related to the band gap temperature-dependence, 
whereas the variation of the high-energy structures orig¬ 
inates from the thermal expansion. 

The results obtained for MgO can be extrapolated 
to other light-element oxides. In the case of XANES 
spectroscopy, the method is applied to corundum in a 
forthcoming publication.®*^ Techniques closely related to 
XANES, such as Non Resonant Inelastic X-ray Scattering 
and core-loss Electron Energy Loss Spectroscopy, may 
also be affected by vibrations and could highly benefit 
from this theoretical framework. 
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Appendix A: Phonon-renormalized X-ray absorption 
cross-section 

In this section we discuss two questions that are rarely 
addressed in the literature. To describe the first ques¬ 
tion recall that, in optical spectroscopy of molecules, 
the equilibrium nuclear positions of the initial and fi¬ 
nal states of the absorption process are different. This 
leads to vibronic effects that are often calculated by using 
the Franck-Condon principle. In X-ray absorption spec¬ 
troscopy, the experimental spectrum is well reproduced 
although the final states are calculated for the same nu¬ 
clear positions as the initial state. How can this be? 

The second question has to do with the effect of tem¬ 
perature on X-ray absorption spectra. This effect is 
often calculated by averaging spectra calculated over 
a distribution of atomic positions corresponding to the 
temperature.^'^^'^® In principle, this procedure is not cor¬ 
rect although it gives reasonable results in practice. How 
can this be? 

In the present section, we answer these two questions 
by showing that the usual calculation methods amount 
to neglecting the nuclear kinetic energy in the absorption 
process and we show how it is possible to go beyond this 
approximation. 

To describe the effect of nuclear vibrations on XANES 
spectra, we write the electric dipole absorption cross- 
section in terms of the wave functions I'l':)^) involving both 
the electronic and the nuclear variables: 

atot (lio;) = ^ I I r 14-)!) I' (5(- A® - fiw) 

= -Airaohu (4'!]| 7^ImG(E° + huj)T , (Al) 
where we used:®^ 

^ 5{Ei - A® - huj) ini = -^ [G(A0 + fitc)] . 

n 

Equation (Al) is equivalent to Eq. (22) for the transition 
operator T = e • r = ' g ' The Green function 

G of the full electronic and nuclear Hamiltonian is the 
solution of the following equation 

[z - Hbo - Tn] G{z) = 1, (A2) 

with the complex energy z = if -|- iy, where 7 an infinites¬ 
imal positive number or a finite number representing the 
broadening due to core-hole lifetime and experimental 
resolution. A straightforward expansion in Eq. (A2) gives 

G(z) = Go(z) + Go{z)TnG{z), (A3) 

where Go(z) = [z — is the Green function in the 

Born-Oppenheimer approximation. Equation (A3) can 
be expanded into 

G(z) = Go(z) + Go{z)TnGo{z) + ... . (A4) 

The resolution of X-ray absorption spectra is determined 
by the lifetime of the core hole (0.36 eV for the A"-edge 
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of Mg®^). It may be considered that the core hole life¬ 
time will smooth out effects that involve a much smaller 
energy. Since the nuclear kinetic energy, evaluated us¬ 
ing the zero-point energy values,®^ is 0.02 eV, it seems 
reasonable to neglect the kinetic energy and to keep 
only the first term of the expansion: 

G{z) ^ Go{z). (A5) 


The second term is expected to play a role in the presence 
of forbidden transitions, but this effect will not be con¬ 
sidered in the present paper. We reach the approximate 
absorption cross-section 

atotihio) = -iiraohuj Im«|r^ Go r|T°). (A6) 


In this approximation, the calculation of the phonon- 
renormalized XANES cross section requires only the 
Born-Oppenheimer Green function Gq, for which we now 
give the following convenient expression: 


(r',R'|Go |r,R) = ^ 


(r^R'|vl/^„) (^jr,R) 

2 ; - e„(R) 


(A7) 


The validity of this expression can be established by 
showing that it solves the equation for Gq- 

[z - Hbo] Go |f,R) = S(r'-K) Sir'-f). (A8) 

Equation (A6) is now evaluated as follows 


atotifizj) = —iiraohuj J drdr'd'Rd'R' Im('I'Q|7~^ I^^R^) 

X (r',R'|Go|f,R)(r,R|r|0. (A9) 

This expression can be simplified by noticing that, in the 
Born-Oppenheimer approximation. Go is diagonal in the 
nuclear variables. Indeed, by writing the initial and final 
states in the Born-Oppenheimer approximation (Eq. 2) 
we obtain 


(f',R'|Go |r,R) = ^ 

n 


2 - £„(R) 


x^x^*(rV„(R)- 


(AlO) 


The completeness relation of Eq. (7) turns Eq. (AlO) into 


(r',R'|5o|r,R) =J(R-R')^ 


V'n(r';R)V'n(i'; R) 


2 - e„(R) 

= <5(R-R') (f'|go(R)|f), (All) 


where go(R) is the electronic Green function for a system 
where the nuclei are fixed at position R. In other words, 
go(R) is the solution of 


(A12) 


where Hbo is evaluated at the nuclear positions R. In¬ 
troducing Eq. (All) in Eq. (A6) the absorption cross- 
section at the first-order in Gq is obtained: 

atotifuz) = -inaohio 

X (r'| 5 o(R)|f)(r,R|r|vk°). (A13) 

The result of Eq. (A13) implies that the XANES calcu¬ 
lation requires only the energy surface - nuclear config¬ 
uration - of the initial state. Hence, Eq. (A13) justihes 
the use of the ground-state crystallographic structure in 
the initial (without core-hole) and final (with core-hole) 
states when calculating the XANES cross-section. The 
total wave functions can be expressed in the BO approx¬ 
imation (Eq. 2) 


drdr'dR Im(«'°|r^ |r',R) 


Gtotify-jj) = —4:7Taohuj Im 

X J dr dr' '^/’o(^^ R) e* • r' 50 (R) e • f R)- (A14) 

Equation (A14) proves that the effect of thermal vibra¬ 
tions on XANES spectra can be obtained by averaging in¬ 
dividual XANES spectra for nuclear positions R weighted 

_ 2 

by the distribution function Xg(R) computed from the 
ground vibrational mode in the ground state. 

However, Eq. (A14) is expressed in a many-body 
framework, whereas AT-edge XANES spectra are usually 
calculated in a single-electron framework. Since this re¬ 
duction is a classical problem, we just give a sketch of 
the derivation. If we rewrite Eq. (AI4) in terms of wave 
functions, we have to deal with matrix elements such as 
('0n|e • f |'0o)i where {tpo) and \tpn) are A^e-body wave 
functions. If we assume that these wave functions can be 
expressed as Slater determinants, the fact that e • r is a 
single-body transition operator implies:®^ 

(V’n|e-r IV'o) = J dr (j)p,ir) e-r (Spir), (A15) 

where i4>p^(j)p') is the only pair of one-electron orbitals 
that are different in {-ipo) and IV'n), where (pp is occupied 
in \ipo) and in |'0n)- For a AT-edge the resulting ex¬ 
pression is 


J dR |xo(R) 


atotifizj) = —47rao^ Im J dR |xo(R)|^ 

X J drdr' 4’ui'r', R) e* • r' ^(r', r; R) e • r (j)isir', R), 

(A16) 

where g(r',r;R) = (r'; R)(/>/ 3 '(r; R )/(2 - e^/). A 

similar expression can be obtained from more sophisti¬ 
cated many-body perturbation theory. 

Considering the cross section in a given nuclear config¬ 
uration from Eq. (21) gives 

o'iot(fiw) = J dR |xo(R)f tT(fc 2 ; R) 


(2 - Hbo) (r'l 5 o(R) |r) = d(r' - r). 


(A17) 





12 


and we demonstrate, restricting ourselves to the first or¬ 
der in the expansion of G(z), that to account for the 
nuclear motion in the XANES cross section one must av¬ 
erage the individual configuration spectra using a proba¬ 
bility distribution, which is consistent with Eq. (19). We 


used a ground-state phonon wave function Xg (R) for no- 
tational convenience. The generalization to a Boltzmann 
distribution p(R) of phonon states at finite temperature 
is straightforward and amounts to replacing |xg(R)| by 
p(R) in Eq. (A17). 
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